
##%%%%%%%%%%%%%%%%%%%%%%%% simOfDemingSIMN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%% Simulation of Deming Data %%%%%%%%%%%%%%%%%%%%%%%
## In this analysis we base our simulations on Deming's data set
## using R package "GenOrd" which allows us to generate samples from 
## correlated ordinal data.  Although the Deming data is actually a 
## mixture of ordinal and continuous data, if we bin the continous 
## variables in sufficiently small groups the difference shouldn't be
## that difference.

##%%%%%%%%%%%%%%% Install latest version of R %%%%%%%%%%%%%%%%%%%%%%%
## One can do this manually by going the R website
## https://cran.r-project.org/bin/windows/base/
## If one has a modern version of R, one can use the "installr" package
## > install.packages("installr") 
## > library(installr) # install+load installr
## Then there will be a button that will be a button on the menu
## [installr] which will allow you to install the latest version of R
## and update all your packages.

##%%%%%%%%%%%%%%%%%%%%%%% Packages %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
## If a package is in the existing library (for example the package
## "foreign") then this can be loaded into the working session using
## > load(foreign)
## If the package isn't in the existing library it can be downloaded
## from a mirror server using
## > install.packages("moments",repos= "http://cran.us.r-project.org")
## > library(moments)
## For this to work, the permissions on the library in the home library
## must be set to allow 'write' for "users".
## Packages can also be loaded using the drop down menu "Packages"
## by clicking on "Packages" button then clicking on the selected
## package name.
## NOTE!! These installations only work if the working directory is
## left as the default.
## One can go back to the basic directory using:
## > setwd(.libPaths()[1])
## If the working directory has been changed to a local directory using:
## > setwd(WorkDir)

library(GenOrd)
library(Matrix)
library(MASS) ## For ginv() - Moore penrose inverse


##%%%%%%%% Set Working Directory %%%%%%%%%%%%
## !!! Here one needs to set one's own working directory
## where the script and data files are held !!!
WorkDir="D:/RWORK/Repository Files"
setwd(WorkDir)

##%%%%%%%%%% Read in function files %%%%%%%%%
source("FUNC.R")

##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%% Start Simulation Section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## Reload the saved files
load("tins2")
load("tins3")
load("tins4")
load("tins4S")
load("tins7")
load("tins7S")

##%%%%%%%%%%%%% Make the data ordinal %%%%%%%%%%%%

## With only a couple of exceptions these variables are hardly
## changed from the originals
tins2O=makeOrdinalM(tins2)
tins3O=makeOrdinalM(tins3)
tins4O=makeOrdinalM(tins4)
tins4SO=makeOrdinalM(tins4S)
tins7O=makeOrdinalM(tins7)
tins7SO=makeOrdinalM(tins7S)


##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%% Simulations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## In this section we will produce simulated data based using the 'ordsample' 
## function from the GenOrd package.
## ordsample(n, marginal, Sigma, support = list(), Spearman = FALSE,
##           cormat = "discrete")


##%%%% Sim for tins2
simInput2 = calcSupportM(tins2O)
simInputF2=supportCorrFix(simInput2)

contSigma2=ordcont(marginal=simInputF2$marginal, Sigma=simInputF2$Sigma,
            support=simInputF2$support,epsilon = 0.37, maxit = 400)[[1]] 
## Didn't work at 0.36, did at 0.37

##%%%% Sim for tins3
simInput3 = calcSupportM(tins3O)
simInputF3=supportCorrFix(simInput3)

contSigma3=ordcont(marginal=simInputF3$marginal, Sigma=simInputF3$Sigma,
            support=simInputF3$support,epsilon = 0.01, maxit = 400)[[1]]

## Didn't work at 0.001, did work at 0.01

##%%%% Sim for tins4
simInput4 = calcSupportM(tins4O)
simInputF4=supportCorrFix(simInput4)

contSigma4=ordcont(marginal=simInputF4$marginal, Sigma=simInputF4$Sigma,
            support=simInputF4$support,epsilon = 0.37, maxit = 400)[[1]]

## Didn't work at 0.35 did work at 0.37

##%%%% Sim for tins4S
simInput4S = calcSupportM(tins4SO)
simInputF4S=supportCorrFix(simInput4S)

contSigma4S=ordcont(marginal=simInputF4S$marginal, Sigma=simInputF4S$Sigma,
            support=simInputF4S$support,epsilon = 0.37, maxit = 400)[[1]]

## Didn't work at 0.35, did work at 0.37

##%%%% Sim for tins7
simInput7 = calcSupportM(tins7O)
simInputF7=supportCorrFix(simInput7)

contSigma7=ordcont(marginal=simInputF7$marginal, Sigma=simInputF7$Sigma,
            support=simInputF7$support,epsilon = 0.28, maxit = 400)[[1]]

## Didn't work at 0.25 did work at 0.28

##%%%% Sim for tins7S
simInput7S = calcSupportM(tins7SO)
simInputF7S=supportCorrFix(simInput7S)

contSigma7S=ordcont(marginal=simInputF7S$marginal, Sigma=simInputF7S$Sigma,
            support=simInputF7S$support,epsilon = 0.27, maxit = 400)[[1]]

## Didn't work at 0.25, did work at 0.27


##%%%%%%%%%%%%%% Save these if you like %%%%%%%%%%%%%%%%%%% 
# save(simInputF2,file="simInputF2")
# save(contSigma2,file="contSigma2")
# save(simInputF3,file="simInputF3")
# save(contSigma3,file="contSigma3")
# save(simInputF4,file="simInputF4")
# save(contSigma4,file="contSigma4")
# save(simInputF4S,file="simInputF4S")
# save(contSigma4S,file="contSigma4S")
# save(simInputF7,file="simInputF7")
# save(contSigma7,file="contSigma7")
# save(simInputF7S,file="simInputF7S")
# save(contSigma7S,file="contSigma7S")

##%%%%%%%%%% Generate large simulation sets %%%%%%%%%%%%%%%%%
simL2 = ordsample(n=1000000, marginal=simInputF2$marginal, Sigma=contSigma2,
                                 support=simInputF2$support,cormat = "continuous")
simL3 = ordsample(n=1000000, marginal=simInputF3$marginal, Sigma=contSigma3,
                                 support=simInputF3$support,cormat = "continuous")
simL4 = ordsample(n=1000000, marginal=simInputF4$marginal, Sigma=contSigma4,
                                 support=simInputF4$support,cormat = "continuous")
simL4S = ordsample(n=1000000, marginal=simInputF4S$marginal, Sigma=contSigma4S,
                                 support=simInputF4S$support,cormat = "continuous")
simL7 = ordsample(n=1000000, marginal=simInputF7$marginal, Sigma=contSigma7,
                                 support=simInputF7$support,cormat = "continuous")
simL7S = ordsample(n=1000000, marginal=simInputF7S$marginal, Sigma=contSigma7S,
                                 support=simInputF7S$support,cormat = "continuous")

## Define targets and inputs for the simulation
t1L2=simL2[,1]
t2L2=simL2[,2]
insL2=simL2[,3:ncol(simL2)]
t1L3=simL3[,1]
t2L3=simL3[,2]
insL3=simL3[,3:ncol(simL3)]
t1L4=simL4[,1]
t2L4=simL4[,2]
insL4=simL4[,3:ncol(simL4)]
t1L4S=simL4S[,1]
t2L4S=simL4S[,2]
insL4S=simL4S[,3:ncol(simL4S)]
t1L7=simL7[,1]
t2L7=simL7[,2]
insL7=simL7[,3:ncol(simL7)]
t1L7S=simL7S[,1]
t2L7S=simL7S[,2]
insL7S=simL7S[,3:ncol(simL7S)]

##%%%%%%%%%%%%%%% Save These %%%%%%%%%%%%%%%%%%%%%%%%
save(t1L2,file="t1L2");save(t2L2,file="t2L2");save(insL2,file="insL2");
save(t1L3,file="t1L3");save(t2L3,file="t2L3");save(insL3,file="insL3");
save(t1L4,file="t1L4");save(t2L4,file="t2L4");save(insL4,file="insL4");
save(t1L4S,file="t1L4S");save(t2L4S,file="t2L4S");save(insL4S,file="insL4S");
save(t1L7,file="t1L7");save(t2L7,file="t2L7");save(insL7,file="insL7");
save(t1L7S,file="t1L7S");save(t2L7S,file="t2L7S");save(insL7S,file="insL7S");



##%%%%% Test the simulation produced models against originals %%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

##%%%%%%%%%%%%%% Build models based on tins data %%%%%%%%%%%%%%%%%%
modt12=lm(tins2[,1]~tins2[,3:ncol(tins2)])
sumModt12=summary(modt12)
modt13=lm(tins3[,1]~tins3[,3:ncol(tins3)])
sumModt13=summary(modt13)
modt14=lm(tins4[,1]~tins4[,3:ncol(tins4)])
sumModt14=summary(modt14)
modt14S=lm(tins4S[,1]~tins4S[,3:ncol(tins4S)])
sumModt14S=summary(modt14S)
modt17S=lm(tins7S[,1]~tins7S[,3:ncol(tins7S)])
sumModt17S=summary(modt17S)
modt17=lm(tins7[,1]~tins7[,3:ncol(tins7)])
sumModt17=summary(modt17)

##%%%%%%%%%%% Compute input sets for simulations %%%%%%%

## for N=1125
Nfe=1125
IIfe = 1:1125
Ivecfe=(0:799)*1125
Imat1=matrix(0,1125,800)
for(i in 1:800){ Imat1[,i]=IIfe+Ivecfe[i]}

## for N=500
N=500
IvecS=(0:1999)*500
IIS = 1:500
Imat.5=matrix(0,500,2000)
for(i in 1:2000){ Imat.5[,i]=IIS+IvecS[i]}

##%%%%%%%%%%%% Do testing %%%%%%%%%%%%%%%%%%
coefMatt12=matrix(0,length(modt12$coef),10)
coefMatt13=matrix(0,length(modt13$coef),10)
coefMatt14=matrix(0,length(modt14$coef),10)
coefMatt14S=matrix(0,length(modt14S$coef),10)
coefMatt17S=matrix(0,length(modt17S$coef),10)
coefMatt17=matrix(0,length(modt17$coef),10)

for(i in 1:10){
	ii=Imat1[,i]
	coefMatt12[,i]=j2z(lm(t1L2[ii]~insL2[ii,])$coef)
	coefMatt13[,i]=j2z(lm(t1L3[ii]~insL3[ii,])$coef)
	coefMatt14[,i]=j2z(lm(t1L4[ii]~insL4[ii,])$coef)
	coefMatt14S[,i]=j2z(lm(t1L4S[ii]~insL4S[ii,])$coef)
	coefMatt17S[,i]=j2z(lm(t1L7S[ii]~insL7S[ii,])$coef)
	coefMatt17[,i]=j2z(lm(t1L7[ii]~insL7[ii,])$coef)
	}

tins2L=cbind(matrix(1,nrow(tins2),1),tins2[,3:ncol(tins2)])
tins3L=cbind(matrix(1,nrow(tins3),1),tins3[,3:ncol(tins3)])
tins4L=cbind(matrix(1,nrow(tins4),1),tins4[,3:ncol(tins4)])
tins4SL=cbind(matrix(1,nrow(tins4S),1),tins4S[,3:ncol(tins4S)])
tins7SL=cbind(matrix(1,nrow(tins7S),1),tins7S[,3:ncol(tins7S)])
tins7L=cbind(matrix(1,nrow(tins7),1),tins7[,3:ncol(tins7)])


fitt12=tins2L%*%coefMatt12
fitt13=tins3L%*%coefMatt13
fitt14=tins4L%*%coefMatt14
fitt14S=tins4SL%*%coefMatt14S
fitt17S=tins7SL%*%coefMatt17S
fitt17=tins7L%*%coefMatt17

corMat=matrix(0,6,10)
corMat[1,]=cor(modt12$fitted.values,fitt12)
corMat[2,]=cor(modt13$fitted.values,fitt13)
corMat[3,]=cor(modt14$fitted.values,fitt14)
corMat[4,]=cor(modt14S$fitted.values,fitt14S)
corMat[5,]=cor(modt17S$fitted.values,fitt17S)
corMat[6,]=cor(modt17$fitted.values,fitt17)

Rmean(corMat)

# > Rmean(corMat)
#            sum
# [1,] 0.7985097
# [2,] 0.9841039
# [3,] 0.8697022
# [4,] 0.9270182
# [5,] 0.8307176
# [6,] 0.7604343

## The best correspondence is found for fit3, fit4, fit4S
## with reliance on SES variables.  Those without, fit2, 
## fit7, fit7S doen't match as well.

##%%%% Try testing correlation to models fit on several
## simulation data sets to another model on a seperate
## simulation data set. 
ii=Imat1[,1]
insL2L=cbind(matrix(1,length(ii),1),insL2[ii,])
insL3L=cbind(matrix(1,length(ii),1),insL3[ii,])
insL4L=cbind(matrix(1,length(ii),1),insL4[ii,])
insL4SL=cbind(matrix(1,length(ii),1),insL4S[ii,])
insL7SL=cbind(matrix(1,length(ii),1),insL7S[ii,])
insL7L=cbind(matrix(1,length(ii),1),insL7[ii,])


fitt12x=insL2L[ii,]%*%coefMatt12
fitt13x=insL3L[ii,]%*%coefMatt13
fitt14x=insL4L[ii,]%*%coefMatt14
fitt14Sx=insL4SL[ii,]%*%coefMatt14S
fitt17Sx=insL7SL[ii,]%*%coefMatt17S
fitt17x=insL7L[ii,]%*%coefMatt17


corMat2=(corMat-corMat)[,2:10]
corMat2[1,]=cor(fitt12x)[1,2:10]
corMat2[2,]=cor(fitt13x)[1,2:10]
corMat2[3,]=cor(fitt14x)[1,2:10]
corMat2[4,]=cor(fitt14Sx)[1,2:10]
corMat2[5,]=cor(fitt17Sx)[1,2:10]
corMat2[6,]=cor(fitt17x)[1,2:10]

# > Rmean(corMat2)
#           sum
# [1,] 0.6177331
# [2,] 0.9793580
# [3,] 0.7798403
# [4,] 0.8840449
# [5,] 0.7354853
# [6,] 0.6600074

## The best correspondence is found for fit3, fit4, fit4S
## again. This shows that the the correspondence pretty much
## matches that seen comparing correlation between the original
## fitted values and simulated coefficient matrices.  Apparently, 
## the lower corresponding values have a higher natural model variance.

## Collect these results in a matrix
corMatA=round(cbind(Rmean(corMat),Rmean(corMat2)),3)

rownames(corMatA)=c("L2","L3","L4","L4S","L7S","L7")
colnames(corMatA)=c("mean corr. orig.fit vs. sim.fit",
				"mean corr. sim.fit1 vs. sim.fit9-10")

# > corMatA
#     mean corr. orig.fit vs. sim.fit mean corr. sim.fit1 vs. sim.fit9-10
# L2                            0.799                               0.618
# L3                            0.984                               0.979
# L4                            0.870                               0.780
# L4S                           0.927                               0.884
# L7S                           0.831                               0.735
# L7                            0.760                               0.660

# write.csv(corMatA,file="simVsOrg.csv")


##%%%%%% fit models for all original data %%%%%%%%
modt12=lm(tins2[,1]~tins2[,3:ncol(tins2)])
sumModt12=summary(modt12)
modt13=lm(tins3[,1]~tins3[,3:ncol(tins3)])
sumModt13=summary(modt13)
modt14=lm(tins4[,1]~tins4[,3:ncol(tins4)])
sumModt14=summary(modt14)
modt14S=lm(tins4S[,1]~tins4S[,3:ncol(tins4S)])
sumModt14S=summary(modt14S)
modt17S=lm(tins7S[,1]~tins7S[,3:ncol(tins7S)])
sumModt17S=summary(modt17S)
modt17=lm(tins7[,1]~tins7[,3:ncol(tins7)])
sumModt17=summary(modt17)

modt22=lm(tins2[,2]~tins2[,3:ncol(tins2)])
sumModt22=summary(modt22)
modt23=lm(tins3[,2]~tins3[,3:ncol(tins3)])
sumModt23=summary(modt23)
modt24=lm(tins4[,2]~tins4[,3:ncol(tins4)])
sumModt24=summary(modt24)
modt24S=lm(tins4S[,2]~tins4S[,3:ncol(tins4S)])
sumModt24S=summary(modt24S)
modt27S=lm(tins7S[,2]~tins7S[,3:ncol(tins7S)])
sumModt27S=summary(modt27S)
modt27=lm(tins7[,2]~tins7[,3:ncol(tins7)])
sumModt27=summary(modt27)


##%%%%%%%% Define II mats for following cv computations
## Here the in and out data are the same 
## !! ignore out/outL.stat's !!
iiin=1:1115
iiout=1:1115
iioutL=1:1115


##%%%%%%%% Compute rsqMat's for all the original models %%%%%%%
rsqMt12=makeRsqMat0.1(tins2[,1],tins2[,3:ncol(tins2)],iiin,iiout,iioutL)
rsqMt22=makeRsqMat0.1(tins2[,2],tins2[,3:ncol(tins2)],iiin,iiout,iioutL)

rsqMt13=makeRsqMat0.1(tins3[,1],tins3[,3:ncol(tins3)],iiin,iiout,iioutL)
rsqMt23=makeRsqMat0.1(tins3[,2],tins3[,3:ncol(tins3)],iiin,iiout,iioutL)


rsqMt14=makeRsqMat0.1(tins4[,1],tins4[,3:ncol(tins4)],iiin,iiout,iioutL)
rsqMt24=makeRsqMat0.1(tins4[,2],tins4[,3:ncol(tins4)],iiin,iiout,iioutL)

rsqMt14S=makeRsqMat0.1(tins4S[,1],tins4S[,3:ncol(tins4S)],iiin,iiout,iioutL)
rsqMt24S=makeRsqMat0.1(tins4S[,2],tins4S[,3:ncol(tins4S)],iiin,iiout,iioutL)

rsqMt17=makeRsqMat0.1(tins7[,1],tins7[,3:ncol(tins7)],iiin,iiout,iioutL)
rsqMt27=makeRsqMat0.1(tins7[,2],tins7[,3:ncol(tins7)],iiin,iiout,iioutL)

rsqMt17S=makeRsqMat0.1(tins7S[,1],tins7S[,3:ncol(tins7S)],iiin,iiout,iioutL)
rsqMt27S=makeRsqMat0.1(tins7S[,2],tins7S[,3:ncol(tins7S)],iiin,iiout,iioutL)

## Merge these
rsqNamesN1=c("rsqMt12","rsqMt22","rsqMt13","rsqMt23","rsqMt14","rsqMt24",
		"rsqMt14S","rsqMt24S","rsqMt17S","rsqMt27S","rsqMt17","rsqMt27")

origNamesN1=c("xx2t1N1","xx2t2N1","xx3t1N1","xx3t2N1","xx4t14N1","xx4t2N1",
		"xx4St1N1","xx4St2N1","xx7St1N1","xx7St2N1","xx7t1N1","xx7t2N1")

## Get column sets
rsqCols=colnames(rsqMt14)[1:59]
rsqColsNoOut=rsqCols[substring(rsqCols,1,3)!="out"]
jjModStats=c(1:6,9:12)

## Merge
foo=get(rsqNamesN1[1])
foo1=foo[1,jjModStats]
rownames(foo)=paste(origNamesN1[1],rownames(foo))
rsqOrigN1=foo[,rsqColsNoOut]
modStatsOrigN1=foo1


for(i in 2:length(rsqNamesN1)){
	foo=get(rsqNamesN1[i])
	foo1=foo[1,jjModStats]
	rownames(foo)=paste(origNamesN1[i],rownames(foo))
	rsqOrigN1=rbind(rsqOrigN1,foo[,rsqColsNoOut])
	modStatsOrigN1=rbind(modStatsOrigN1,foo1)
	}
rownames(modStatsOrigN1)=origNamesN1

## Test to see if "t value" and "tstatF" match
# > max(abs(abs(rsqOrigN1[,"t value"])-rsqOrigN1[,"tstatF"]))
# [1] 2.002711e-07
##%%%%%%%%%%%%%%% rsqMat's for "Male's" %%%%%%%%%%%%%%%%%%
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
malE=iiin[tins2[,"Male"]==1]

malEin=malE
malEout=malE
malEoutL=malE

rsqMt12m=makeRsqMat0.1(tins2[,1],tins2[,3:ncol(tins2)],malEin,malEout,malEoutL)
rsqMt22m=makeRsqMat0.1(tins2[,2],tins2[,3:ncol(tins2)],malEin,malEout,malEoutL)

rsqMt13m=makeRsqMat0.1(tins3[,1],tins3[,3:ncol(tins3)],malEin,malEout,malEoutL)
rsqMt23m=makeRsqMat0.1(tins3[,2],tins3[,3:ncol(tins3)],malEin,malEout,malEoutL)


rsqMt14m=makeRsqMat0.1(tins4[,1],tins4[,3:ncol(tins4)],malEin,malEout,malEoutL)
rsqMt24m=makeRsqMat0.1(tins4[,2],tins4[,3:ncol(tins4)],malEin,malEout,malEoutL)

rsqMt14Sm=makeRsqMat0.1(tins4S[,1],tins4S[,3:ncol(tins4S)],malEin,malEout,malEoutL)
rsqMt24Sm=makeRsqMat0.1(tins4S[,2],tins4S[,3:ncol(tins4S)],malEin,malEout,malEoutL)

rsqMt17m=makeRsqMat0.1(tins7[,1],tins7[,3:ncol(tins7)],malEin,malEout,malEoutL)
rsqMt27m=makeRsqMat0.1(tins7[,2],tins7[,3:ncol(tins7)],malEin,malEout,malEoutL)

rsqMt17Sm=makeRsqMat0.1(tins7S[,1],tins7S[,3:ncol(tins7S)],malEin,malEout,malEoutL)
rsqMt27Sm=makeRsqMat0.1(tins7S[,2],tins7S[,3:ncol(tins7S)],malEin,malEout,malEoutL)



## Merge this data
rsqNamesN.5=c("rsqMt12m","rsqMt22m","rsqMt13m","rsqMt23m","rsqMt14m","rsqMt24m",
		"rsqMt14Sm","rsqMt24Sm","rsqMt17Sm","rsqMt27Sm","rsqMt17m","rsqMt27m")

origNamesN.5=c("xx2t1N.5","xx2t2N.5","xx3t1N.5","xx3t2N.5","xx4t1N.5","xx4t2N.5",
		"xx4St1N.5","xxt4St2N.5","xx7St1N.5","xx7St2N.5","xx7t1N.5","xx7t2N.5")


## Get column sets
rsqCols=colnames(rsqMt14)[1:59]
rsqColsNoOut=rsqCols[substring(rsqCols,1,3)!="out"]
jjModStats=c(1:6,9:12)

## Merge
foo=get(rsqNamesN.5[1])
foo1=foo[1,jjModStats]
rownames(foo)=paste(origNamesN.5[1],rownames(foo))
rsqOrigN.5=foo[,rsqColsNoOut]
modStatsOrigN.5=foo1


for(i in 2:length(rsqNamesN.5)){
	foo=get(rsqNamesN.5[i])
	foo1=foo[1,jjModStats]
	rownames(foo)=paste(origNamesN.5[i],rownames(foo))
	rsqOrigN.5=rbind(rsqOrigN.5,foo[,rsqColsNoOut])
	modStatsOrigN.5=rbind(modStatsOrigN.5,foo1)
	}
rownames(modStatsOrigN.5)=origNamesN.5

## Test to see if "t value" and "tstatF" match
# > max(abs(abs(rsqOrigN.5[,"t value"])-rsqOrigN.5[,"tstatF"]))
# [1] 3.657965e-09



##%%%%%%%%%%% Reduced tables for paper %%%%%%%%%%%%
foo=modStatsOrigN1[,"cv.r.squaredI"]/modStatsOrigN1[,"r.squaredI"]
modStatsOrigN1s=cbind(modStatsOrigN1[,c(1:3,5,7)],foo)
colnames(modStatsOrigN1s)[ncol(modStatsOrigN1s)]="cv.rsq/in.rsq"

## Put all t1's and t2's together
cc=c(3:1,4:6)
modStatsOrigN1s=modStatsOrigN1s[,cc]

# > round(modStatsOrigN1s[1:6,],3)
#          n-inputs n-obs numb. inputs w. pval<=0.1 r.squaredI cv.r.squaredI cv.rsq/in.rsq
# xx2t1N1        62  1115                        12      0.156         0.048         0.306 rsqQ2
# xx3t1N1        10  1115                         5      0.174         0.158         0.907 rsqQ5
# xx4t14N1       69  1115                        15      0.242         0.130         0.537 rsqQ3
# xx4St1N1       35  1115                        12      0.213         0.157         0.735 rsqQ4
# xx7St1N1       28  1115                        10      0.071         0.019         0.266 rsqQ2
# xx7t1N1        62  1115                        20      0.126         0.011         0.084 rsqQ1


# > round(modStatsOrigN1s[7:12,],3)
#        n-inputs n-obs numb. inputs w. pval<=0.1 r.squaredI cv.r.squaredI cv.rsq/in.rsq
# xx2t2N1        62  1115                        11      0.114         0.008         0.069 rsqQ1
# xx3t2N1        10  1115                         7      0.092         0.075         0.812 rsqQ4
# xx4t2N1        69  1115                        11      0.153         0.042         0.272 rsqQ2
# xx4St2N1       35  1115                         9      0.115         0.055         0.477 rsqQ2
# xx7St2N1       28  1115                         5      0.069         0.015         0.213 rsqQ2
# xx7t2N1        62  1115                        11      0.134         0.017         0.128 rsqQ1

## The same thing for N.5
foo=modStatsOrigN.5[,"cv.r.squaredI"]/modStatsOrigN.5[,"r.squaredI"]
modStatsOrigN.5s=cbind(modStatsOrigN.5[,c(1:3,5,7)],foo)
colnames(modStatsOrigN.5s)[ncol(modStatsOrigN.5s)]="cv.rsq/in.rsq"

## Put all t1's and t2's together
## reorder columns
cc=c(3:1,4:6)
modStatsOrigN.5s=modStatsOrigN.5s[,cc]

# > round(modStatsOrigN.5s[1:6,],3)
#           n-inputs n-obs numb. inputs w. pval<=0.1 r.squaredI cv.r.squaredI cv.rsq/in.rsq
# xx2t1N.5        56   557                         9      0.213         0.008         0.036 rsqQ1
# xx3t1N.5         9   557                         2      0.165         0.134         0.815 rsqQ4
# xx4t1N.5        63   557                        12      0.283         0.078         0.274 rsqQ2
# xx4St1N.5       34   557                         9      0.232         0.114         0.494 rsqQ2
# xx7St1N.5       28   557                         4      0.072        -0.042        -0.582 rsqQ1
# xx7t1N.5        62   557                         8      0.136        -0.179        -1.311 rsqQ1

# > round(modStatsOrigN.5s[7:12,],3)
#            n-inputs n-obs numb. inputs w. pval<=0.1 r.squaredI cv.r.squaredI cv.rsq/in.rsq
# xx2t2N.5         56   557                         5      0.124        -0.113        -0.907 rsqQ1
# xx3t2N.5          9   557                         3      0.087         0.055         0.634 rsqQ3
# xx4t2N.5         63   557                         7      0.179        -0.067        -0.373 rsqQ1
# xxt4St2N.5       34   557                         4      0.131         0.002         0.013 rsqQ1
# xx7St2N.5        28   557                         4      0.090        -0.023        -0.259 rsqQ1
# xx7t2N.5         62   557                         8      0.162        -0.093        -0.575 rsqQ1

## Save these for section on original data
write.csv(modStatsOrigN1s,file="modStatsOrigN1s.csv")
write.csv(modStatsOrigN.5s,file="modStatsOrigN.5s.csv")





